La régression géographiquement pondérée : GWR
Comment prendre en compte l’effet local du spatial en statistique
Cette fiche présente la réalisation d’une analyse de données à l’aide de la régression géographique pondérée ou GWR. La modélisation statistique “classique” a beaucoup de mal à gérer le spatial. Et pour cause la dimension spatiale va à l’encontre d’une règle fondamentale de la statistique : l’indépendance. Or lors d’une analyse d’un phénomène social il est souvent mauvais de considérer la dimension spatiale d’une données comme un simple aléa dans le meilleur des cas ou pire de l’ignorer complétement. L’objectif de cette fiche est de vous présenter une méthode qui vous permettra d’étudier concrétement l’effet du spatial dans un modèle de régression. C’est à dire de mesurer statiquemet l’effet de l’information spatiale sur une variable à expliquer. Cette analyse sera reproductible, les données spatiales utilisées proviennent de la base ADMIN-EPXRESS-COG de l’IGN. Les données statitiques elles ont été construites à partir de la base des notaires de France, il s’agit de données de recherches mises à disposition.
-
Prérequis : connaissance de base en analyse de données
statistiques, avec au moins une compréhension de ce qu’est une
corrélation et une régression?.
Pourquoi la GWR
En statistique la méthode reine pour étudier, analyser la nature des relations et des effets entre des variables est le modèle de régression linéaire. Le principe de la régression linéaire est de modéliser la variable que nous souhaitons étudier (aussi appeler variable dépendante, VD) comme une fonction linéaire des variables que nous aurons définis comme explicatives de la VD (aussi appelées variables indépendantes, VI). Cette fonction linéaire nous permets d’obtenir des coefficient (appelés beta) et des résidus (noté epsilon). Ces betas représentent l’effet de nos VI sur notre VD. Ces betas sont considérés comme globaux, autrement dit sans variation.
Or lorsque l’on s’intéresse à un phénomène social avec une emprise sur un espace, un territoire donnée cette méthode pose deux problèmes majeurs :
1- Le premier est empirique, Les recherches en sciences sociales et notamment en géographie montrent qu’il est très rare qu’un phénomène social soit constant dans l’espace. C’est le concept d’hétérogénéité spatiale, l’effet de nos VI va varier en fonction de l’espace. Ainsi, un coefficient qui soit globale et uniforme pour mesurer un effet est tentant mais pas tenable, sur ce point nous pouvons nous référer à l’artcile de Brundson et al. (1996). Ce concept d’hétérogénéité dans l’espace ce traduit en statistique par celui de non stationarité.
2- Le deuxième est statistique : les tests statistiques doivent répondre à un certain nombre de conditions et particulièrement une régression linéaire. Ainsi, pour être valide certain critère doivent être validé. On pense notamment à la notion d’indépendance et d’autocorrélation des résidus. Or de part leur nature les données spatiales ne peuvent pas remplir ces conditions fondamentale pour une régression classique. C’est la 1er loi de la géographie de Tobler : “everything is related to everything else, but near things are more related than distant things”
La GWR va nous permettre ainsi de résoudre ces deux problèmes en intégrant la dimension spatiale de nos données tout en tenant compte de l’hétérogénéité (ou non staionarité) de leur effet.
Les packages
Voici les packages que nous utiliserons :
# Chargement, visualisation et manipulation de la données
library(here)
library(DT)
library(dplyr)
# Analyse et réprésentation statistique
library(car)
library(correlation)
library(corrplot)
library(ggplot2)
library(gtsummary)
library(GGally)
library(plotly)
# Manipulation et représentation de la données spatiales (cartographie)
library(sf)
library(mapsf)
library(rgeoda)
library(RColorBrewer)
# Calcul du voisinage et réalisation de la GWR
library(spdep)
library(GWmodel)Vous pouvez vérifier l’installation des différents packages à l’aide des lignes de codes suivantes.
# Packages nécessaires
my_packages <- c("here", "DT", "dplyr", "car", "correlation", "corrplot", "ggplot2", "gtsummary", "GGally",
"plotly", "sf", "mapsf", "rgeoda", "RColorBrewer", "spdep", "GWmodel")
# Vérifier si ces packages sont installés
missing_packages <- my_packages[!(my_packages %in% installed.packages()[,"Package"])]
# Installation des packages manquants depuis le CRAN
if(length(missing_packages)) install.packages(missing_packages,
repos = "http://cran.us.r-project.org")
# Chargement des packages nécessaires
lapply(my_packages, library, character.only = TRUE)1 Présentation et préparation des données
Pour réaliser cette fiche nous aurons à la fois besoin de données statistiques et de données spatiales.
Comme indiqué en introduction les données spatiales proviennent de la base ADMIN-EPXRESS-COG de l’IGN, l’unité sera ici l’EPCI et le format initial le shapefile.
Nos donnéees statistiques seront sous le format csv. Il s’agit de données construites par Frédéric Audard et Alice Ferrari depuis la base de données des notaires de France. Ce fichier a été simplifié pour ne conserver que les variables d’intérêts parmi une cinquantaine.
1.1 Chargement des données sur le prix de l’immobilier par EPCI
library(here)
# On situe le dossier dans lequel se trouve nos données
csv_path <- here("data", "donnees_standr.csv")
immo_df <- read.csv2(csv_path)
#Pour visualiser les données dans le doc
datatable(head(immo_df, 10))Ce fichier est composé des 10 variables suivantes :
- SIREN : code SIREN de l’EPCI
- prix_med : pris médian par EPCI à la vente (au m2 ?)
- perc_log_vac : % logements vacants
- perc_maison : % maisons
- perc_tiny_log : % petits logements (surface < ?)
- dens_pop : densité de population (nb habitants / km2 ?)
- med_niveau_vis : médiane du niveau de vie
- part_log_suroccup : % logements suroccupés
- part_agri_nb_emploi : % agriculteurs
- part_cadre_profintellec_nbemploi : % cadres et professions intellectuelles
La variable SIREN nous servira de “clés” pour joindre
ces données statistiques aux données spatiale, la variable
prix_med sera la variable que nous chercherons à expliquer
(VD) toutes les autres seront nos variables explicatives (VI)
Lorsque l’on s’apprête à faire de la modélisation statistique il est très recommandé de réaliser cette opération au moins sur les variables que vous utiliserez comme variables explicatives dans votre modèle.
Sur R on peut facilement réaliser cette opération avec la fonction scale(). Ou à “la main” l’opération est simple. On soustrait chaque valeur par la moyenne puis on divise par l’écart-type.
1.2 Chargement des données géographiques : les EPCI de France métropolitaine
Ces données proviennent de la base ADMIN-EPXRESS-COG de l’IGN, édition 2021. Le format d’entrée est le shapefile mais nous passerons par une conversion au format sf, ce qui nous permet d’utiliser le package mapsf, pour les prévisualiser :
library(here)
library(sf)
shp_path <- here("data", "EPCI.shp")
epci_sf <- st_read(shp_path)Reading layer `EPCI' from data source
`/mnt/9A8C6DB98C6D9115/gregoire/GWR_rzine/data/EPCI.shp' using driver `ESRI Shapefile'
Simple feature collection with 1242 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 124110 ymin: 6049642 xmax: 1242428 ymax: 7110453
Projected CRS: RGF93_Lambert_93
library(mapsf)
# pour voir les données géographiques
mf_map(x = epci_sf)# et la table attributaire correspondante
datatable(head(epci_sf, 10))1.3 Jointure des données géographiques et tabulaires
Les 2 données n’ont pas le même nombre de lignes :
nrow(immo_df)[1] 1223
nrow(epci_sf)[1] 1242
On constate nos deux jeux de données n’ont pas exactement le même
nombre de ligne. En effet, le jeux de données immo_df
possède moins de ligne que notre oject sf epci_sf. Cela
indique simplement que nous n’avons pas l’indication du prix médian de
l’immobilier pour tous les epci de France métropolitaine. Il peut être
intéressant d’identifier et visualiser les EPCI qui n’ont pas de
correspondance dans le tableau de données immo_df, pour ce
faire on réalise la jointure on conservant toutes les lignes de
epci_sf :
# l'option all.x = TRUE permet de garder toutes les lignes de epci_sf,
# même celles qui n'ont pas de correspondance dans immo_df
data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN", all.x = TRUE)
nrow(data_immo)[1] 1242
# on peut filtrer les données de la jointure pour ne voir que les epci n'ayant pas de correspondance dans le tableau immo_df
datatable(data_immo[is.na(data_immo$prix_med),])mf_map(x = data_immo[is.na(data_immo$prix_med),])Cependant, notre VD étant prix_med les lignes vides ne
nous intéressent pas, nous ne les conserverons pas car pourraient poser
problèmes lors de la réalisation de nos analyses. On refait la jointure
en ne gardant que les EPCI ayant une correspondance dans le tableau de
données :
data_immo <- merge(x = epci_sf, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN")
nrow(data_immo)[1] 1223
datatable(head(data_immo, 10))Dans le cas où vous préféreriez manipuler vos données sous un format shape, ou dans le cas où ce format serait requis pour utiliser certains packages ou certaines formules, vous pouvez convertir votre objet sf en shape à l’aide de la ligne de commande suivante :
shape <- as(data_immo, "Spatial")2 Création du voisinage
Avant de procéder à nos différentes analyse nous devons d’abord créer et définir notre voisinage. Cette étape est absolument essentielle. En effet, cette notion de voisinage est centrale en statistique spatiale, le principe de base étant que le voisinage a un effet sur nos individus. Les choix qui seront fait dans la construction du voisinage impacteront de fait très fortement les résultats.
2.1 Voisinage
Nous ne développerons pas ici tout ce qu’est et ce qu’implique la définition d’un voisinage. Pour cela, nous vous renvoyons vers les travaux de XXXX ou l’ouvrage XXXXX.
Ce qui est important ici de savoir c’est qu’un voisinage peut être de 3 type :
- Basé sur la contigüté,
- Basé sur la distance
- Basé sur la proximité
Lorsque nous travaillons avec des polygones (comme c’est le cas ici) le plus souvent on va se baser sur une matrice de contiguité. Il faut encore savoir qu’il existe plusieurs type de voisinage basé sur la contigüté. Dans un cas classique nous utiliserons celui de type queen. Queen est une référence à la reine des échecs. La reine au échecs peut se déplacer dans toutes les directions, et ici on va considérer les voisins contigus à notre polygone de tous côtés. Il s’oppose au type rook qui fait référence à la tour, les voisins seront donc définis à partir des mouvements de cette pièce sur l’échiquier.
Figure 2.8 du manuel INSEE Codifier la structure de voisinage (insee?)
Heureusement R permet assez simplement de définir notre voisinage.
library(spdep)
# Création de la liste des voisins : avec l'option queen = TRUE,
# sont considérés comme voisins 2 polygones possédant au moins 1 sommet commun
#help(poly2nb)
neighbours_epci <- poly2nb(data_immo, queen = TRUE)
# Obtention des coordonnées des centroïdes
coord <- st_coordinates(st_centroid(data_immo))
# cette opération se fait aussi avec un shape
#shape_nbq <- poly2nb(shape, queen = TRUE)
#coord<-coordinates(shape)Voici la réprésentation graphique de notre voisinage.
# Faire un graphe de voisinage
plot(neighbours_epci, coord)Pour comprendre ce que contient neighbours_epci :
# si on prend le 1er élément de neighbours_epci, on voit qu'il a pour voisins les poygones 62, 74 etc.
neighbours_epci[[1]][1] 62 74 338 1135 1136 1137 1140
# ce qu'on peut vérifier sur la carte :
neighbours1 <- data_immo[c(1,62,74,338,1135,1136,1137,1140),]
neighbours1$index <- rownames(neighbours1)
mf_map(x = neighbours1)
mf_label(x = neighbours1, var = "index")Nous précisons qu’un voisinage peut aussi tout à fait se calculer
lorsque l’on a pas de polygones mais simplement des coordonnées (des
points). Les matrices de distance sont alors souvent plus adaptées. Pour
définir le voisinage il faut utiliser les fonctions
knearneigh() et knn2nb()
2.2 Création de la matrice de voisinage
Une fois le voisnage définit nous pouvons créer une matrice de voisinage, qui permettra d’attribuer un poids à chaque voisins.
# la fonction nb2listw attribue des poids à chaque voisin
# par ex. si un polygone a 4 voisins, ils auront chacun un poids de 1/4 = 0.25
#help("nb2listw")
neighbours_epci_w <- nb2listw(neighbours_epci)
# les poids sont stockés dans le 3ème élément de neighbours_epci_w
# par ex. si on veut connaître les poids des voisins du 1er élément :
neighbours_epci_w[[3]][1][[1]]
[1] 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571 0.1428571
# cet élément a 7 voisins qui ont donc un poids de 1/7 soit ~0.14Comment faire si l’on a un polygone sans voisins ? Sur le plan
technique, la fonction nb2listw prévoit ce cas de figure.
Il faut utiliser l’argument zero.policy est indiquer
TRUE.
Au niveau théorique, c’est moins clair. De
manière générale les indices d’autocorrélation spatiale et autres
régression spatiale ont été conçus en partant du principe que les
entités spatiales avaient un voisinage. Ceci dit il n’y a pas à ma
connaissance de régles absolues qui obligent à les supprimer ou
intégrer.
3 Approche statistique “classique”
Nous pouvons donc commencer notre analyse. Avant de se diriger vers une GWR, il est important de suivre la voie “classique”. Elle permet d’abords de mieux connaitre et appréhender nos données mais surtout vérifier que la méthode statitisuqe classique ne parvient pas à rendre compte de la complexité du phénomène sans tenir compte de la dimension spatiale et donc justifier l’usage de la GWR.
Traditionnellement cela passe par trois étape :
1- Etude de la distribution des données. 2- Analyse des corrélations 3- Réalisation d’un odèle linéaire classique
N’étant pas le sujet de la fiche nous passerons rapidement sur ces étapes et nous ne attarderons pas sur des considérations théoriques. Toutefois nous indiquerons des manières de les réaliser sur R et des éléments de lectures.
3.1 Exploration des variables
Cette première étape très importante permet d’étudier la distribution de nos données et d’identifier d’éventuels individus extrêmes qui pourrait venir perturber les résultats.
Ici par exemple nous serions en droit de nous poser la question de conserver dans notre jeux de données l’entité spatiale avec un prix médian à plus de 10000 qui se détache très clairement de tous les autres EPCI.
Pour visualiser la distribution d’une variable quantitative
l’histogramme est une bonne solution. Pour le réaliser nous utilisons
dans ce cas le package plotly qui permet l’intéractivité de
la figure.
# Distribution de la variable dépendante :
library(plotly)
add_histogram(plot_ly(data_immo, x = ~prix_med))Ce que nous faisons avec notre VD il est important de le faire aussi pour nos VI.
# Distribution des variables indépendantes :
a <- add_histogram(plot_ly(data_immo, x = ~log(perc_log_vac), name = "perc_log_vac"))
b <- add_histogram(plot_ly(data_immo, x = ~log(perc_maison), name = "perc_maison"))
c <- add_histogram(plot_ly(data_immo, x = ~log(perc_tiny_log), name = "perc_tiny_log"))
d <- add_histogram(plot_ly(data_immo, x = ~log(dens_pop), name = "dens_pop"))
e <- add_histogram(plot_ly(data_immo, x = ~log(med_niveau_vis), name = "med_niveau_vis"))
f <- add_histogram(plot_ly(data_immo, x = ~log(part_log_suroccup), name = "part_log_suroccup"))
g <- add_histogram(plot_ly(data_immo, x = ~log(part_agri_nb_emploi), name = "part_agri_nb_emploi"))
h <- add_histogram(plot_ly(data_immo, x = ~log(part_cadre_profintellec_nbemploi), name = "part_cadre_profintellec_nbemploi"))
fig = subplot(a, b, c, d, e, f, g, h, nrows = 2)
fig3.2 Etude des corrélations
Pour un point plus complet sur les correlations et leurs mise en oeuvre sur R, vous pouvez vous rendre sur Rzine et la fiche Analyse des corrélations avec easystats, Guide pratique avec R.
Très rapidement, la corrélation permet d’étudier le lien, la relation entre deux variables. La corrélation repose sur la covariance entre les variables. Quand 2 variables covarient, un écart à la moyenne d’une variable est accompagné par un écart dans le même sens (corrélation positive) ou dans le sens opposé de l’autre (corrélation négative) pour le même individus de l’autre variable. Dis autrement, lorsque deux variables covarient pour chaque valeur qui s’écarte de la moyenne, on s’attend à trouver un écart à la moyenne pour l’autre variable.
La conception d’un modèle statistique doit absolument être le fruit d’une réflexion portant sur le choix des variables indépendantes (explicatives) et le choix de la méthode de régression. Et pour définir un modèle de régression certaine régle doivent être respectées.
L’étude des corrélations peut donc apporter une aide précieuse
dans cette réflexion. Elle pourra nous aider dans le choix des variables
à intégrer au modéle mais dans le même temps de vérifier certaines des
conditions de réalisation de notre régression.
Ainsi, une analyse des corrélation pourra vérifier :
l’existence d’un lien entre les variables indépendantes et la variable à étudier. En effet dans une régression linéaire, il est nécessaire d’avoir une relation linéaire entre la VD et les différentes VI.
la multicolinéarité des variables indépendantes. Les corrélations ne doivent pas être trop fortes entre nos VI. Un coefficient > 0.7 en valeurs absolues doit entrainer la suppression des variables concernés. Cela peut aussi être vérifié très efficacement avec le VIF (Variance Inflation Factor) mais peut se faire seulement après avoir lancé notre modèle.
L’absence des corrélation entre les variables explicatives du modèle et les variables externes. En effet, les variables d’influence doivent être incluses dans le modèle (sauf dans le cas où cela induirait une trop grande multicolinéarité).
Pour calculer une matrice de corrélation :
library(correlation)
data_cor <- immo_df %>% select(-SIREN)
immo_cor <- correlation(data_cor, redundant = TRUE)
summary(immo_cor)# Correlation Matrix (pearson-method)
Parameter | part_cadre_profintellec_nbemploi | part_agri_nb_emploi | part_log_suroccup | med_niveau_vis | dens_pop | perc_tiny_log | perc_maison | perc_log_vac | prix_med
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
prix_med | 0.55*** | -0.45*** | 0.60*** | 0.59*** | 0.52*** | 0.49*** | -0.60*** | -0.68*** | 1.00***
perc_log_vac | -0.31*** | 0.39*** | -0.28*** | -0.46*** | -0.21*** | -0.19*** | 0.35*** | 1.00*** |
perc_maison | -0.58*** | 0.54*** | -0.79*** | -0.24*** | -0.46*** | -0.75*** | 1.00*** | |
perc_tiny_log | 0.75*** | -0.51*** | 0.87*** | 0.22*** | 0.62*** | 1.00*** | | |
dens_pop | 0.59*** | -0.29*** | 0.62*** | 0.18*** | 1.00*** | | | |
med_niveau_vis | 0.43*** | -0.33*** | 0.15*** | 1.00*** | | | | |
part_log_suroccup | 0.65*** | -0.43*** | 1.00*** | | | | | |
part_agri_nb_emploi | -0.57*** | 1.00*** | | | | | | |
part_cadre_profintellec_nbemploi | 1.00*** | | | | | | | |
p-value adjustment method: Holm (1979)
# On peut aussi visualiser les corrélation à l'aide d'un corrélogramme
# 1- Méthode simple
# immo_cor %>%
# summary(redundant = FALSE) %>%
# plot(type="tile", show_labels =TRUE, show_p = TRUE, digits = 1, size_text=3)
# 2- Méthode plus complexe mais meilleure visualisation
mat_cor_comp <- summary(immo_cor, redundant = TRUE)
# Nom des lignes = valeurs de la première colonne ("Parameter")
rownames(mat_cor_comp ) <- mat_cor_comp[,1]
# Transformation du data.frame en objet matrice (+ suppression première colonne)
mat_cor<- as.matrix(mat_cor_comp[,-1])
# Calcul du nombre total d'individus
nb <- nrow(data_cor)
# Calcul des matrices de p-values et des intervalles de confiance
p.value <- cor_to_p(mat_cor, n = nb, method = "auto")
# Extraction de la matrice des p-value uniquement
p_mat <- p.value$p
library(corrplot)
library(RColorBrewer)
corrplot(mat_cor,
p.mat = p_mat,
type = "upper",
order = "hclust",
addCoef.col = "white",
tl.col = "gray",
number.cex = 0.5,
tl.cex= 1,
tl.srt = 45,
col=brewer.pal(n = 8, name = "PRGn"),
sig.level = 0.000001,
insig = "blank",
diag = FALSE, )L’études des corrélation nous permet ici de confirmer une relation entre notre variable à expliquer et toutes les variables explicatives définies. Elle met aussi à jour de très fortes multicolinéarités, ce qui va nous poser problème. Dans le cadre de la définition d’un modéle linéaire classique il faudrait sortir du modèle les variables explicatives qui corrèle trop fortement. Dans le cadre de cette fiche nous faisons le choix de conserver toutes nos VI, cela fera sens au moment de la GWR.
3.3 Régression linéaire ou Méthode des moindre carrés ordinaire (MCO).
La régression linéaire est un des modèles les plus utilisés en SHS, elle peut être simple (une seule variable explicative) ou multiple (plusieurs VI). Le principe n’est en réalité pas si complexe. La régression linéaire consiste à modeliser la covariation entre une variable à expliquer et une ou des variables explicatives.
Pour ce faire le modèle de régression va chercher à estimer les termes de l’equation de la droite de régression entre la variable à expliquer et la variable explicative. Cette equation va prendre la forme \(f(x)= ax + b + e_i\). Equation qui ressemble à la fonction affine \(f(x)= ax + b\), qui est étudié en général au collège.
Si l’on cherche à modéliser la covariation entre le prix médian et le pourçentage de logement vacant l’équation de notre droite de régression ressemblerait à : \(prixmed_i = a*perclogvac_i + b + e_i\).
Où \(a\) est le coefficient associé à la variable du pourcentage de logements vacants, \(b\) est la constante qui apparaitra sous le nom d’intercept dans les résultats et enfin \(e\) qui lui correspond aux résidus. Les résidus étant ce qui incarne l’écart au modèle.
Traditionnelement on va faire usage d’une régression linéaire lorsque l’on veut prédire les valeurs de notre VD dans les cas où elle n’aurait pas été mésuré. Ou si l’on souhaite comprendre les relations statistiques entre les variables.
Pour lancer notre modèle de régression linéaire avec toutes nos VI voici le slignes de commandes :
# Dans le fonctionnement sur R il est important de sotcker la régression dans un objet.
# Pour lancer la régression on va utiliser la fonction lm() dont les 2 lettres sont l'acronyme pour linear model
mod.lm <- lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi,
data = data_immo)
# On affiche les principaux résultats avec la fonction summary
summary(mod.lm)
Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log +
dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi +
part_cadre_profintellec_nbemploi, data = data_immo)
Residuals:
Min 1Q Median 3Q Max
-1566.8 -220.2 -27.2 174.0 3277.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1592.873 11.204 142.171 < 2e-16 ***
perc_log_vac -287.337 14.134 -20.330 < 2e-16 ***
perc_maison -128.255 20.041 -6.400 2.22e-10 ***
perc_tiny_log -261.556 27.934 -9.363 < 2e-16 ***
dens_pop 173.942 15.272 11.389 < 2e-16 ***
med_niveau_vis 288.017 13.860 20.781 < 2e-16 ***
part_log_suroccup 388.982 27.352 14.221 < 2e-16 ***
part_agri_nb_emploi -20.785 15.059 -1.380 0.168
part_cadre_profintellec_nbemploi -7.841 19.904 -0.394 0.694
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 391.8 on 1214 degrees of freedom
Multiple R-squared: 0.7784, Adjusted R-squared: 0.7769
F-statistic: 532.9 on 8 and 1214 DF, p-value: < 2.2e-16
#Pour visualiser les résultats de manière plus agréable on peut aussi utiliser la fonction tbl_regression du package gtsummary.
library(gtsummary)
mod.lm %>%
tbl_regression(intercept = TRUE)| Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | 1,593 | 1,571, 1,615 | <0.001 |
| perc_log_vac | -287 | -315, -260 | <0.001 |
| perc_maison | -128 | -168, -89 | <0.001 |
| perc_tiny_log | -262 | -316, -207 | <0.001 |
| dens_pop | 174 | 144, 204 | <0.001 |
| med_niveau_vis | 288 | 261, 315 | <0.001 |
| part_log_suroccup | 389 | 335, 443 | <0.001 |
| part_agri_nb_emploi | -21 | -50, 8.8 | 0.2 |
| part_cadre_profintellec_nbemploi | -7.8 | -47, 31 | 0.7 |
|
1
CI = Confidence Interval
|
|||
# On peut également visualiser graphiquement les coefficients des variables explicatives
GGally::ggcoef_model(mod.lm)3.3.1 Interpréter les résultats
Pour interpréter les résultats plusieurs éléments fournient par la
commande summary(mod.lm) sont importants.
D’abord l’information fournit par Adjusted R-squared, il
s’agit du \(R^2\) qui est le
coefficient de détermination. Il est ici de 0.77 ce qui veut dire que
77% de la variation du prix médian est expliqué par notre modèle. Juste
en dessous, on a le p-value associé à notre modèle. S’il est inférieur à
\(0.05\) on peut considérer qu’il est
statistiquement significatif. Voici pour les infos associé globalement
au modéle.
Ensuite, on peut regarder ce qu’il se passe au niveau des variables
explicatives. La première colonne appelé estimateest le
coefficient de régression associé à la variable explicative. Le signe va
être très important car va donner la direction de la relation
(exactement comme pour les corrélations).Ici pour le pourcentage de
logements vacants il est de \(-287\) ce
qui veut dire que lorsque ce pourcentage augment de \(1\%\) alors le prix median baisse de \(287€\). La colonne \(Pr(>|t|)\) correspond à la p-value
associé à ce résultat. Une fois encore s’il est inférieur à \(0.05\) alors on peut considérer comme
statistiquement significatif. Dans notre cas on peut dire que la part
d’agriculteurs, de cadre et profession intellectuel dans le nombre
d’emploi n’ont pas un effet significatif sur le prix médian du
logement.
La \(t-value\) est le coefficient divisé par son erreur standard. Plus ce quotient est proche de \(0\) plus on peut considérer qu’il n’y a pas d’effet de notre VI sur notre VD. En revanche sans le cas où on aurait plus de 200 individus si $t-value > |1.96| $ alors on peut considérer qu’il y a \(95\%\) de chance qu’il y ait un effet significatif de la VI sur la VD.
Ainsi, le pourcentage de logements vacant à un effet sur le prix médian toutes choses égales quant aux autres variables explicatives.
Pour aller au delà de ce raisonnement il faut intégrer des effets d’intéraction au modèle.
3.4 Diagnostic de notre modèle linéaire
3.4.1 Multicolinéarité
Un des enjeux les plus important dans le cadre de régression multiple est de vérifier la multicolinéarité entre les variables explicatives. Le risque d’une trop grande colinéarité est de biaiser le modèle et notamment de biaiser les estimations des erreurs type des coefficient de régression (et donc les t-value et p-value).
La VIF (Variance Inflation Factor) est une très bonne méthode pour vérifier les risques de multicolinéarité. Elle suppose simplement d’avoir estimer un premier modèle pour être utilisée.
library(car)
vif(mod.lm) perc_log_vac perc_maison
1.577984 3.204058
perc_tiny_log dens_pop
6.221631 1.864236
med_niveau_vis part_log_suroccup
1.532403 5.977951
part_agri_nb_emploi part_cadre_profintellec_nbemploi
1.812419 3.162576
# On peut aussi directement l'ajouter au résumé des coefficient obtenu avec gtsummary
library(gtsummary)
mod.lm %>%
tbl_regression(intercept = TRUE) %>% add_vif()| Characteristic | Beta | 95% CI1 | p-value | VIF1 |
|---|---|---|---|---|
| (Intercept) | 1,593 | 1,571, 1,615 | <0.001 | |
| perc_log_vac | -287 | -315, -260 | <0.001 | 1.6 |
| perc_maison | -128 | -168, -89 | <0.001 | 3.2 |
| perc_tiny_log | -262 | -316, -207 | <0.001 | 6.2 |
| dens_pop | 174 | 144, 204 | <0.001 | 1.9 |
| med_niveau_vis | 288 | 261, 315 | <0.001 | 1.5 |
| part_log_suroccup | 389 | 335, 443 | <0.001 | 6.0 |
| part_agri_nb_emploi | -21 | -50, 8.8 | 0.2 | 1.8 |
| part_cadre_profintellec_nbemploi | -7.8 | -47, 31 | 0.7 | 3.2 |
|
1
CI = Confidence Interval, VIF = Variance Inflation Factor
|
||||
Sur le seuil de VIF à ne pas dépasser les sources en SHS varient en fonction des disciplines, certaines proposent 5 et d’autres même 10. Malgré tout, en géographie le consensus est autour d’une valeur critique de 4. Un VIF supérieur à 4 devrait entrainer la supprésion de la variable du modèle car implique une forte colinéarité et donc un risque élevé de biaiser le modèle. A partir de 3 il convient de s’interroger sur la présence de la variable.
Au delà de la suppression ou non des variables concernées, il est aussi très important de pouvoir identifier les variables à VIF élevés.
On peut facilement représenter graphiquement les scores de VIF.
library(car)
score_vif <- vif(mod.lm)
# On peut aussi directement l'ajouter au résumé des coefficient obtenu avec gtsummary
barplot(score_vif, main = "VIF Values", horiz = TRUE, col = "steelblue", las=2)
#ajout du seuil de 4
abline(v = 4, lwd = 3, lty = 2)
# et de la limite de 3
abline(v = 3, lwd = 3, lty = 2)Comme le laisser supposer l’étude des corrélations de nos variables, nous avons en effet une forte multicolinéarité entre certaines de nos variables explicatives. Selon le VIF nous devrions donc relancer notre modèle sans les variables fortement colinéaires, c’est à dire sans le pourcentage de logements sur occupé et sans le pourcentage de petits logements. Pour commencer, on peut retirer du modèle la variable ayant le VIF le plus élevé à savoir le pourcentage de petits logements.
mod.lm2 <- lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop +
med_niveau_vis + part_log_suroccup + part_agri_nb_emploi +
part_cadre_profintellec_nbemploi, data = data_immo)
summary(mod.lm2)
Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop +
med_niveau_vis + part_log_suroccup + part_agri_nb_emploi +
part_cadre_profintellec_nbemploi, data = data_immo)
Residuals:
Min 1Q Median 3Q Max
-1502.2 -226.3 -42.3 173.5 3535.9
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1592.549 11.597 137.329 < 2e-16 ***
perc_log_vac -328.301 13.911 -23.601 < 2e-16 ***
perc_maison -100.010 20.507 -4.877 1.22e-06 ***
dens_pop 161.783 15.750 10.272 < 2e-16 ***
med_niveau_vis 280.586 14.322 19.591 < 2e-16 ***
part_log_suroccup 237.069 22.792 10.401 < 2e-16 ***
part_agri_nb_emploi 2.696 15.370 0.175 0.861
part_cadre_profintellec_nbemploi -77.744 19.098 -4.071 4.99e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 405.5 on 1215 degrees of freedom
Multiple R-squared: 0.7623, Adjusted R-squared: 0.761
F-statistic: 556.8 on 7 and 1215 DF, p-value: < 2.2e-16
vif(mod.lm2) perc_log_vac perc_maison
1.426783 3.131469
dens_pop med_niveau_vis
1.850756 1.527378
part_log_suroccup part_agri_nb_emploi
3.874580 1.762159
part_cadre_profintellec_nbemploi
2.717634
library(gtsummary)
mod.lm2 %>%
tbl_regression(intercept = TRUE) %>% add_vif()| Characteristic | Beta | 95% CI1 | p-value | VIF1 |
|---|---|---|---|---|
| (Intercept) | 1,593 | 1,570, 1,615 | <0.001 | |
| perc_log_vac | -328 | -356, -301 | <0.001 | 1.4 |
| perc_maison | -100 | -140, -60 | <0.001 | 3.1 |
| dens_pop | 162 | 131, 193 | <0.001 | 1.9 |
| med_niveau_vis | 281 | 252, 309 | <0.001 | 1.5 |
| part_log_suroccup | 237 | 192, 282 | <0.001 | 3.9 |
| part_agri_nb_emploi | 2.7 | -27, 33 | 0.9 | 1.8 |
| part_cadre_profintellec_nbemploi | -78 | -115, -40 | <0.001 | 2.7 |
|
1
CI = Confidence Interval, VIF = Variance Inflation Factor
|
||||
GGally::ggcoef_model(mod.lm2)On note qu’au niveau global il y a peu de changement le \(R^2\) a très légérement baissé, on passe d’une variation expliqué à 77% à un taux d’explication de 76% et le modèle est toujours significatif. Les changement les plus important se situent au niveau des effets partiels. L’effet de la part de cadre et de profession intellectuel dans le nombre d’emploi sur le prix médian du logement est devenu significatif et on constate même que le VIF du pourcentage de logement sur occupé est passé sous le seuil critique.
3.4.2 Principe de Parcimonie
Lorsque l’on conçoit un modèle linéaire, nous sommes sencés respecter un principe de parcimonie. Ce principe implique qu’un bon modèle à un nombre optimal de variable. Bref, qu’il ne s’embarasse pas de variables non significatives. Ce principe veut donc que nous retirons de notre modèle la variable part d’agriculteurs dans le nombre d’emploi.
La fonction step() permet de réaliser une régression pas
à pas descendante (ou ascendante). Dans le cas d’une regression
descendante, le modèle initial comprend toutes les variables, comme pour
la régression standard mais cette fois l’algorithme va retirer la
variable ayant la plus faible contribution au modèle si la variation du
\(R^2\) n’est pas significative en
l’éliminant. La procédure va être répétée jusqu’à ce que toutes les
variables conservées contribuent significativement à l’amélioration du
R2. La régression descendante va donc retirer les variables non
significatives une à une. Ainsi, le dernier modèle proposé doit contenir
toutes les variables ayant une contribution significative au \(R^2\).
# L'argument "both" permeet d'utiliser les deux méthodes : ascendante et ascendante
step(mod.lm2, direction = "both")Start: AIC=14696.68
prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis +
part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi
Df Sum of Sq RSS AIC
- part_agri_nb_emploi 1 5060 199817567 14695
<none> 199812507 14697
- part_cadre_profintellec_nbemploi 1 2725360 202537867 14711
- perc_maison 1 3911245 203723752 14718
- dens_pop 1 17351111 217163618 14796
- part_log_suroccup 1 17792190 217604698 14799
- med_niveau_vis 1 63121115 262933622 15030
- perc_log_vac 1 91601521 291414028 15156
Step: AIC=14694.71
prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis +
part_log_suroccup + part_cadre_profintellec_nbemploi
Df Sum of Sq RSS AIC
<none> 199817567 14695
+ part_agri_nb_emploi 1 5060 199812507 14697
- part_cadre_profintellec_nbemploi 1 3211704 203029272 14712
- perc_maison 1 4155413 203972980 14718
- dens_pop 1 17567092 217384659 14796
- part_log_suroccup 1 18007058 217824626 14798
- med_niveau_vis 1 63117884 262935451 15028
- perc_log_vac 1 94962190 294779757 15168
Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop +
med_niveau_vis + part_log_suroccup + part_cadre_profintellec_nbemploi,
data = data_immo)
Coefficients:
(Intercept) perc_log_vac
1592.55 -327.82
perc_maison dens_pop
-99.01 162.05
med_niveau_vis part_log_suroccup
280.55 237.44
part_cadre_profintellec_nbemploi
-78.93
On observe ici qu’à la fin la variable
part_agri_nb_emploin’est donc pas conservé.
Notre nouveau modèle devrait donc ressembler à ça :
mod.lm3 <- lm(formula = prix_med ~ perc_log_vac + perc_maison + dens_pop + med_niveau_vis + part_log_suroccup + part_cadre_profintellec_nbemploi, data = data_immo)3.4.3 Analyser les résidus :
L’analyse des résidus est très importante car les conditions de validité d’un modèle linéaire au delà des résultats repose grandement sur les résidus. Ils permettent en outre d’identifier les individus extrêmes (ou outliers).
Pour rappel, les résidus correspondent à l’écart au modèle. Ainsi, un residu > 0 implique que notre individu a été sous-estimé par le modèle (il est au dessus de la droite de régression), un résidu < 0 que l’individu a été sur-estimé par le modèle (il est sous la droite de régression).
Les 3 conditions qui concernent les résidus sont :
- Ils doivent suivre une loi normale.
- Ils ne doivent pas varier en fonction des variables explicatives. C’est l’hypothèse d’homoscédasticité, ils ont une variance homogène.
- Ils ne doivent pas être autocorrélés.
Soyons clair lorsque la démarche est de simplement réaliser une étude de modèle linéaire il est rare de voir des articles en géographie où ces trois conditions sont étudiés ou validé. C’est pourtant important même s’il faut reconnaitre que les types de données en SHS remplissent pas toujours ces conditions. Ceci dit, dans une démarche qui s’arreterait au modèle linéaire s’il y en a une à vérifier ce serait plutôt la normalité des résidus.
Pour obtenir les résidus :
res_modlm <- mod.lm$residuals
datatable(as.data.frame(res_modlm))On peut également les visualiser :
par(mfrow=c(1,3))
qqPlot(mod.lm) #diagramme quantile-quantile qui permet de vérifier l'ajustement d'une distribution à un modèle théorique, ici loi normal[1] 36 266
hist(rstudent(mod.lm), breaks = 50, col="darkblue", border="white", main="Analyse visuelle des résidus") # Histogramme pour donner une autre indication sur la normalité
plot(rstudent(mod.lm)) # un graphique pour visualiser l'homoscédasticité des résidusSi la voie graphique ne vous inspire pas il existe des tests statistique qui permettent de vérifier la normalité des résidus ou bien leur homoscédasticité. Ils ont cela de particulier qu’ici nous cherchons à accepter H0 et donc pour valider la normalité ou l’homoscédasticité il faut que \(p-value>0.05\)
# Pour étudier la normalité on peut utiliser le test de Shapiro-Wilk
shapiro.test(mod.lm$residuals)
Shapiro-Wilk normality test
data: mod.lm$residuals
W = 0.90792, p-value < 2.2e-16
# Pour évaluer l'homoscdasticité on peut utiliser le test de Breusch-Pagan. Le package car propose une fonction pour le réaliser
ncvTest(mod.lm)Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 1151.844, Df = 1, p = < 2.22e-16
Dans les deux cas il nous faut rejeter H0, les résidus n’ont donc pas une distribution normale et il y a hétéroscdasticité de la variance des résidus.
Le modèle ne peut donc pas être analysé en l’état. Le problème de l’hétéroscédasticité des résidus indique un problème de spécification du modèle (par exemple une variable manquante).
Apparté sur les outliers :
Le qqplot nous indique deux individus extrême ici les ceux ayant pour identifiant 266 et 36. Il peut dans certain cas être intéressant de supprimer ces individus et voir comment réagit le modèle.
# Pour visualiser les individus concernés
data_immo[c(36, 266),]Simple feature collection with 2 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 637302.2 ymin: 6532738 xmax: 1013257 ymax: 6879128
Projected CRS: RGF93_Lambert_93
CODE_SIREN ID NOM
36 200023372 EPCI____0000002150235984 CC de la Vallée de Chamonix-Mont-Blanc
266 200054781 EPCI____0000002150236221 Métropole du Grand Paris
NATURE prix_med perc_log_vac perc_maison perc_tiny_log
36 Communauté de communes 6400 -2.3174835 -2.558919 -0.07390879
266 Métropole 10920 -0.7044846 -3.732690 6.43751449
dens_pop med_niveau_vis part_log_suroccup part_agri_nb_emploi
36 -0.2437935 0.7430146 1.873292 -1.027161
266 22.7709258 0.9077121 7.345980 -1.068415
part_cadre_profintellec_nbemploi geometry
36 0.07587777 MULTIPOLYGON (((997329.6 65...
266 6.30681026 MULTIPOLYGON (((666798.7 68...
# Pour relancer un nouveau modèle sans l'individus le plus extrême. Notez que l'on peut en supprimer plusieurs d'un coup avec subset =-c(36,266)
mod.lmx <- update(mod.lm, subset=-266)
# Etudier le nouveau modèle
summary(mod.lmx)
Call:
lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log +
dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi +
part_cadre_profintellec_nbemploi, data = data_immo, subset = -266)
Residuals:
Min 1Q Median 3Q Max
-1622.6 -215.0 -37.5 170.0 2767.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1587.246 10.526 150.797 < 2e-16 ***
perc_log_vac -302.182 13.317 -22.692 < 2e-16 ***
perc_maison -132.159 18.814 -7.024 3.57e-12 ***
perc_tiny_log -227.353 26.356 -8.626 < 2e-16 ***
dens_pop -4.725 19.978 -0.237 0.813
med_niveau_vis 284.947 13.012 21.899 < 2e-16 ***
part_log_suroccup 404.159 25.701 15.725 < 2e-16 ***
part_agri_nb_emploi -17.545 14.138 -1.241 0.215
part_cadre_profintellec_nbemploi 23.384 18.841 1.241 0.215
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 367.8 on 1213 degrees of freedom
Multiple R-squared: 0.7823, Adjusted R-squared: 0.7809
F-statistic: 544.9 on 8 and 1213 DF, p-value: < 2.2e-16
vif(mod.lmx) perc_log_vac perc_maison
1.589305 3.168447
perc_tiny_log dens_pop
6.073006 2.089258
med_niveau_vis part_log_suroccup
1.531885 5.726830
part_agri_nb_emploi part_cadre_profintellec_nbemploi
1.811307 3.111740
par(mfrow=c(1,2))
plot(rstudent(mod.lmx))
qqPlot(mod.lmx)[1] 36 180
# Il est possible de comparer les deux modèles et les coefficients
car::compareCoefs(mod.lm, mod.lmx, pvals = TRUE)Calls:
1: lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log +
dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi +
part_cadre_profintellec_nbemploi, data = data_immo)
2: lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log +
dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi +
part_cadre_profintellec_nbemploi, data = data_immo, subset = -266)
Model 1 Model 2
(Intercept) 1592.9 1587.2
SE 11.2 10.5
Pr(>|z|) < 2e-16 < 2e-16
perc_log_vac -287.3 -302.2
SE 14.1 13.3
Pr(>|z|) < 2e-16 < 2e-16
perc_maison -128.3 -132.2
SE 20.0 18.8
Pr(>|z|) 1.6e-10 2.1e-12
perc_tiny_log -261.6 -227.4
SE 27.9 26.4
Pr(>|z|) < 2e-16 < 2e-16
dens_pop 173.94 -4.72
SE 15.27 19.98
Pr(>|z|) < 2e-16 0.81
med_niveau_vis 288.0 284.9
SE 13.9 13.0
Pr(>|z|) < 2e-16 < 2e-16
part_log_suroccup 389.0 404.2
SE 27.4 25.7
Pr(>|z|) < 2e-16 < 2e-16
part_agri_nb_emploi -20.8 -17.5
SE 15.1 14.1
Pr(>|z|) 0.17 0.21
part_cadre_profintellec_nbemploi -7.84 23.38
SE 19.90 18.84
Pr(>|z|) 0.69 0.21
Dans certain cas les différences sont mineures, ici la différence est importante, en effet, on voit qu’une VI a perdu sa significativité. Quoi qu’uil en soit c’est une opération qui doit être effectuée avec prudence, la suppression des individus pose toujours question notamment en terme de justification théorique. Il faut le faire uniquement si l’analyse des individus indique un problème important (valeurs aberrantes, inversion…).
3.4.4 L’autocorrélation des résidus
C’est la condition la plus difficile à vérifier et celle qui pose le plus problèmes. Heureusement la géographie c’est doté d’outil pour mesurer notamment l’autocorrélation spatial. En réalité ici nous espérons très fortement qu’il y ait une autocorrélation spatial. Cela rendrait notre modèle linéaire classique caduc mais nous permettrai de justifier l’utilisation de la régression spatiale.
Les deux outils connus au moins de nom par tous les géographes sont les tests de Moran et celui de Geary.
Dans la littérature le test de Moran semble être préféré à celui de Geary en raison d’une stabilité plus grande.
library(spdep)
# test de moran des residus de la régression H0: pas d'autocorrélation spatiale
lm.morantest(model = mod.lm,
listw = neighbours_epci_w)
Global Moran I for regression residuals
data:
model: lm(formula = prix_med ~ perc_log_vac + perc_maison +
perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup +
part_agri_nb_emploi + part_cadre_profintellec_nbemploi, data =
data_immo)
weights: neighbours_epci_w
Moran I statistic standard deviate = 21, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Observed Moran I Expectation Variance
0.359707408 -0.003517610 0.000299161
# Test de Geary H0 pas d'autocorrélation.
# Attention : Pour avoir le coefficient il faut faire 1-"Résultat test de Geary" (soit ici le coefficient est 0.67)
# Le coefficient de Geary s'étend de 0 à 2, 1 étant le "0" et signifiant aucune corrélation. Par ailleurs, un score inférieurs à 1 imlique une corrélation positive et un score supérieur à 1 une corrélation négative.
geary(x = data_immo$prix_med,
listw = neighbours_epci_w,
n = length(neighbours_epci),
n1 = length(neighbours_epci)-1,
S0 = Szero(neighbours_epci_w))$C
[1] 0.3324317
$K
[1] 18.24222
On voit dans les deux cas qu’il y aurait bien une auto-corrélation spatiale. Cela implique deux choses très importantes :
- La condition d’absence d’autocorrélation de nos résidus n’est pas vérifié le modèle classique n’est pas interprétable en l’état.
- La dimension spatiale joue un rôle, nous pouvons justifier d’étudier de manière plus approfondi l’autocorrélation spatiale et faire l’usage de la GWR.
3.4.4.1 Cartographie des residus de la régression
On intégre les résidus à la table attributaire de notre objet sf. A priori, comme on a utilisé nos données spatiale (sf) pour la régression les données sont classées dans le bon ordre.
data_immo$res_reg <- mod.lm$residualsLa carte des résidus :
# Définition d'une palette de couleur
cols_v1 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272", "#fb6a4a", "#de2d26")
# Réalisation de la carte
mf_map(x = data_immo,
var = "res_reg",
type = "choro",
border = "#ebebeb",
lwd = 0.1,
breaks=quantile(data_immo$res_reg,seq(0,1, by=1/11)),
pal=cols_v1,
leg_title = "Résidus de régression\nlinéaire 'classique'",
leg_val_rnd = 1)
mf_title("Résidus modèle lm") #titreSur cette carte on voit très clairment une spatialisation des résidus, sans même faire les tests nous aurions pu voir que la dimension spatiale jouait bien un rôle. Sans autocorrelation nous aurions eu une repartition aléatoire des résidus.
4 Analyse de l’autocorrélation spatiale
La question de l’autocorrélation est centrale, c’est en effet elle qui rend notre modèle linéaire inopérant.
A ce stade, nous pouvons nous demander ce qu’est l’autocorrélation spatiale ? Il s’agit tout simplement de la corrélation, positive ou négative, d’une variable avec elle même du fait de la localisation spatiale des observations. Si l’autocorrélation spatiale est positive mes données seront semblables à celles de mes voisins et dissemblables de celles des individus éloignés. A l’inverse si l’autocorrélation spatiale est négative mes données seront différentes de celles de mes voisins et ressembleront davantages à celles des individus éloignés.
L’étude de l’autocorrélation spatiale est particulièrement intéressante car permet de mieux comprendre l’éventuelle structure spatiale du phénomène observé. C’est d’autant plus important que lorsque qu’il y a effectivement une structure spatiale sous jacente on ne peut généralement pas faire appel à la pluspart des statistiques classiques qui repose sur l’hypothèse d’indépendance, qui du fait de cette dimension spatiale n’est plus respecté.
L’analyse de l’autocorrélation spatiale se fait à deux niveaux :
- Le niveau global
- Le niveau local
4.1 Niveau global
Pour mesure l’autocorrélation comme nous l’avons vue précédemment les deux outils les plus utilisés sont les test de Moran et de Geary. L’indice de Moran va considérer les variances et covariances en prenant compte de la différence entre chaque et la moyenne de toutes les observations. L’indice de Geary de son côté prend en compte la différence entre les observations voisines.
Pourquoi l’un plutôt que l’autre ? Il semblerait que Moran soit jugé plus stable et revienne plus souvent dans les articles scientifiques. Ceci dit le coût de réalisation n’est pas très élevé rien n’empêche de faire l’un et l’autre pour voir si les résultats sont cohérent entre eux.
Commençons par représenter sur une carte notre variable du prix médian des logements par EPCI
# La palette de couleur:
cols_v1 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272", "#fb6a4a", "#de2d26")
# Carte du prix médian
mf_map(x = data_immo,
var = "prix_med",
type = "choro",
border = "#ebebeb",
lwd = 0.1,
breaks=quantile(data_immo$prix_med,seq(0,1, by=1/11)),
pal=cols_v1,
leg_title = "Prix médian",
leg_val_rnd = 1)
mf_title("Prix Médian du logement par EPCI France Métropolitaine") #titre
A l’oeil, on voit assez nettement qu’une structure spatiale semble
exister.
Vérifions le statistiquement !
library(spdep)
#Pour l'occasion on va standardiser notre prix median. Cela permettra par la suite de le comparer à d'autres variables si d'autres analyses d'autocorrélation spatiale sont réalisés
data_immo$prix_med_z<- (data_immo$prix_med-mean(data_immo$prix_med))/sd(data_immo$prix_med)
# Test de Geary
# Attention à la lecture particulière des résultats de l'indice de Geary
geary.test(data_immo$prix_med_z, neighbours_epci_w, zero.policy = TRUE, randomisation = FALSE)
Geary C test under normality
data: data_immo$prix_med_z
weights: neighbours_epci_w
Geary C statistic standard deviate = 35.972, p-value < 2.2e-16
alternative hypothesis: Expectation greater than statistic
sample estimates:
Geary C statistic Expectation Variance
0.3324317242 1.0000000000 0.0003444045
# Test de Moran :
# On indique dans un premier temps la variable que l'on souhaite analyser
# Puis la matrice de voisinnage
# L'argument zero.policy=TRUE permet de préciser que l'on souhaite intégrer à l'analyse les entités spatiales qui n'auraient pas de voisins
# L'argument randomisation=FALSE transmet l'instruction à la fonction que nous supposons que la distribution est normale. Dans le cas contraire on devrait partir sur une solution de type Monte-Carlo qui repose sur la randomisation
moran.test(data_immo$prix_med_z, neighbours_epci_w, zero.policy = TRUE, randomisation = FALSE)
Moran I test under normality
data: data_immo$prix_med_z
weights: neighbours_epci_w
Moran I statistic standard deviate = 41.143, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.7154340217 -0.0008183306 0.0003030652
# Ce test d'autocorrélation se lit exactement comme un test de corrélation classique. On va donc s'interesser au signe et à la grandeur du coefficient et à la p-value du test.
# Ici on donc une autocorrélation spatiale positive et importante.Qu’est ce qu’implique l’existence de cette autocorrélation ? Comme il l’ a été mentioné on pourrait définir l’autocorrélation spatiale comme la corrélation, positive ou négative, d’une variable avec elle-même du fait de la localisation spatiale des observations. Cette autocorrélation spatiale peut, donc être d’une part, être le résultat d’un phenomène inobservé (un aléa) ou qu’on ne peut mesurer qui intervient dans l’espace et qui de fait crée une structure spatiale. Il existe différent phénomène sociaux de la sorte comme par exemple des phénomènes d’interaction ou de diffusion (comme les phénomènes de diffusion technologique). Ces phénomènes peuvent produire de l’autocorrélaton spatiale. D’autre part, dans la définition modèle statistique, la mesure de l’autocorrélation spatiale peut être envisagée comme un outil de diagnostic et de détection d’un “mauvais” (du point de vue statistique) modèle (variables non intégrées spatialement corrélées, erreurs sur le choix de l’échelle à laquelle le phénomène spatial est analysé, etc.)
En bref, le calcul de l’autocorrélation vous permettra soit d’identifier et mettre au jour un phénomène spatial non mésuré ou de vérifier la qualité et iabilité de votre modèle
Pour visualiser rapidement la structure spatiale on peut aussi réaliser un diagramme de Moran qui est complémentaire au test statistique.
moran.plot(data_immo$prix_med_z,neighbours_epci_w, labels = TRUE, xlab="prix medians par epci" , ylab="moyenne du prix médians par epci des voisins")Comment lire ce diagramme ?
Si nos individus sont répartis complétement aléatoirement dans l’espace alors c’est le signe d’une absence d’autocorrélation, la pente de la droite de régression sera donc de 0. Si on contraire la pente est non nulle, comme c’est le cas ici, c’est bien le signe de la présence d’une autocorrélation spatiale.
Un aspect important de ce diagramme c’est qu’il permet d’ores et déjà de caractériser nos coefficients d’autocorrélation spatiale. Comme ce graphique est centré sur la moyenne qui est de 0, tous les points à droite de l’axe des ordonnées auront une moyenne > 0 et ceux à gauche <0. Cette réflexion s’applique également à l’axe des abscisses. Par convention on désigne les individus avec une moyenne > 0 par le terme high et ceux avec une moyenne < 0 par le terme low, au sens de supérieur ou inférieur à la moyenne.
Ainsi on peut découper ce diagramme en 4 quadrants. Les quadrants en haut à droite et en-bas à gauche correspondent aux individus ayant une autocorrélation spatiale positive. C’est à dire une valeur proche de celle de leurs voisins. Pour le quadrant en haut à droite on parle du quadrant High-High composé d’individus ayant une valeur de la variable plus élevée que la moyenne entouré d’autres individus qui lui ressemble. Pour le quadrant en bas à gauche on parle du quadrant Low-Low composé d’individus avec un score plus faible que la moyenne avec des voisins avec un score similaire.
Le quadrant en bas à droite est considéré comme le quadrant High-Low. On y retrouve des individus avec un score plus élevé que la moyenne sur la variable du prix médian mais avec un voisinage qui ne lui ressemble pas. Autocorrélation spatiale négative mais score élevé à la variable. En haut à gauche on retrouve à l’inverse les individus avec une valeur du prix médian plus faible que la moyenne et un indice d’autocorrélation spatiale négatif. C’est le quadrant Low-High.
4.2 Niveau local
Les indice de Geary et Moran ont une limite assez importante. Ils partent du principe que le processus spatial étudié serait stationnaire, autrement dit global. Cela veux dire que l’effet de la dimension serait le même dans tout notre espace. Ce qui pose sérieusement question et ceux d’autant plus que l’emprise géographique augmente. La compréhension et la réalisation d’autocorrélation au niveau locale permet d’avncer par la suite vers la GWR.
C’est Luc Anselin qui va dévellopé ce concept et définir ce concept d’indicateur d’autocorrélation spatiale locale. Il parlera de LISA (Local Indicators of Spatial Association). Luc Anselin définit que le LISA de chaque individus statistique indique l’intensité du regroupement spatial de valeurs similaires autour de cette individu. Dis légérement autrement un inividu avec un LISA élevé va avoir une concentration autour de lui de voisins avec des valeurs similaires( pour nous par exemple un regroupement d’individus avec un prix particulièrement élevé ou à l’inverse particulièrement bas).
Le LISA le plus utilisé est le I de Moran Local.
L’idée de faire appel au LISA et de compléter le niveau global par une approche locale. On va chercher à la fois à détecter des clusters qui correspondent à un regroupement significatif de valeurs identiques dans une zone particulière et de repérer des zones de non stationnarité c’est à dire qui ne suivent pas le processus global.
Le logiciel GéoDa a été développé par Luc Anselin et son équipe notamment pour étudier l’autocorrélation spatiale et les LISA. C’est un très bon logiciel en clic-boutons, avec une documentation riche.
Calcul de LISA sur R avec le I de Moran Local
# Pour ce faire on va utiliser le package rgeoda développé également par Luc Anselin pour réaliser sur R les traitement de GeoDa
library(rgeoda)
# Pour utiliser la fonction local_moran proposé par le package rgeoda 2 pré-requis:
# 1- utiliser la fonction queen_weights du package rgeoda pour calculer une matrice de contiguité de type queen
queen_w <- queen_weights(data_immo)
# 2- Sortir la variable à étudier dans un vecteur
prix_med_z = data_immo["prix_med_z"]
lisa <- local_moran(queen_w, prix_med_z)
# Pour visualiser les résultats des LISA il faut les stocker dans des objets ou dans les base de données pour les représenter
lisa_colors <- lisa_colors(lisa)
lisa_labels <- lisa_labels(lisa)
lisa_clusters <- lisa_clusters(lisa)
lisa_value <- lisa_values(lisa)
lisa_pvalue<- lisa_pvalues(lisa)Pour illustrer les Moran locaux on peut réaliser une carte
## Carte de Moran locaux
data_immo$moranlocalvalue<- lisa_values(lisa)
cols_v2 <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#eff3ff", "#ffffce", "#fee5d9", "#fcbba1", "#fc9272", "#fb6a4a", "#de2d26")
mf_map(x = data_immo,
var = "moranlocalvalue",
type = "choro",
border = "#ebebeb",
lwd = 0.1,
pal=cols_v2,
leg_title = "Local Moran",
leg_val_rnd = 1)
mf_title("Carte des LISA du prix médian du logement")Si le score du Moran local est > 0 cela implique que l’on avoir un regroupement de valeurs similaires, plus faibles ou plus élevées que la moyenne. En revanche si le score est < 0 alors on a un regroupement de valeurs dissimilaires, par exemple des valeurs plus faibles entourés de valeur plus fortes (centre de la Corse).
L’étude des p-value associés est également importante car des LISA significatifs (p-value<0.05) renvoient à des clusters de valeurs (similaires ou dissimilaires) qui serait plus marqué que ce l’on peut observer si la répartition spatiale du phénomène était aléatoire.
# Carte des p-value des moran locaux
data_immo$moranlocalpvalue<- lisa_pvalues(lisa)
# Pour plus de lisibilité dans la carte on va faire des classes des p-value
data_immo<- data_immo %>% mutate(lisapvalue_fac = case_when(moranlocalpvalue <= 0.002 ~ "[0.001;0.002[",
moranlocalpvalue <= 0.01 ~ "[0.002;0.01[",
moranlocalpvalue <= 0.05 ~ "[0.01;0.05[",
TRUE ~ "[0.05;0.5]")) %>%
mutate(lisapvalue_fac = factor(lisapvalue_fac,
levels = c("[0.001;0.002[", "[0.002;0.01[",
"[0.01;0.05[",
"[0.05;0.5]")))
mypal <- mf_get_pal(n = 4, palette = "Reds")
mf_map(x = data_immo,
var = "lisapvalue_fac",
type = "typo",
border = "grey3",
lwd = 0.1,
pal=mypal,
leg_title = "P-value Local Moran")
mf_title("Carte des significativité des LISA")Les regroupements que l’on observe vont pouvoir se rapprocher des 4 type de regroupement que nous avions sur le diagramme de Moran. Cette carte des clusters est représenté selon une convention particulière.
# En utilisant le package mapsf
data_immo$testmoran <- sapply(lisa_clusters, function(x){return(lisa_labels[[x+1]])})
colors <- c("white","blue",rgb(0,0,1,alpha=0.4),rgb(1,0,0,alpha=0.4),"red")
mf_map(x = data_immo,
var = "testmoran",
type = "typo",
border = "black",
lwd = 0.1,
pal= colors,
val_order = c("Not significant","Low-Low","Low-High","High-Low","High-High"),
leg_title = "Lisa cluster")
mf_title("LISA clusters")Si on est en présence d’une autocorrélation spatiale au niveau global, les LISA pourront aussi servir d’indicateur d’instabilité locale. Ils indiquent à la fois des clusters locaux qui vont avoir un impact fort sur le processus spatial global (un score d’autocorrélation locale plus marquée que l’autocorrélation globale) ou à l’inverse des zones que se démarquent du processus global (plus faible autocorrélation).
En revanche, s’il n’y a pas d’autocorrélation spatiale au niveau globale les LISA vont nous permettre de détecter des zones où des valeurs semblables se regroupent de façon significative. Ils font émerger des structures au niveau local au sein desquels les liens entre voisins seront particulièrement important.
5 Régression géographiquement pondérée (GWR)
Bien que l’analyse d’autocorrélation spatiale soit riche en “enseignements” sur nos données il reste cependant une étape qui est celle de la modélisation. Un des intérêts de la régression spatiale au même titre qu’une régression classique va être de mieux comprendre la relation qui unit les variables explicatives à la variable étudiée. En effet, si avec le LISA on a pu éventuellement mettre en évidence un effet local du spatial à l’échelle des EPCI de France métropolitaine nous avons pas encore étudié les effets (et leur variabilité) de nos VI sur notre VD.
Il existe différents modèle de régression spatiale, toute la question est de savoir quel modèle utilisé ? Ce choix va dépendre de la nature des phénomènes spatiaux étudiés.
Pour bien comprendre, un peu de théorie est nécessaire. Luc Anselin va distinguer d’un côté ce qui relève de l’autocorrélation spatiale. Il y a autocorrélation spatiale positive lorsque qu’il y a une similarité entre les valeurs d’une même zone spatiale. Dans ce cas de figure on utilisera des modèle SEM, SDEM, SAR… De l’autre côté il y a l’hétérogénéité, qui renvoie à une instabilité, on aurait une variabilité spatiale de nos paramètres. L’idée c’est que nos VI peuvent avoir un effet qui n’est pas le même partout dans l’espace. Dans ce cas nous opterons pour la régression géographiquement pondérée (GWR).
Ainsi, si on s’intéresse au lien entre les voisins on est dans l’autocorrélation spatiale et les modèle SEM, SAR & cie mais si on s’intéresse à l’hétérogénéité de nos variable c’est à dire à leur variabilité selon leur localisation on est dans la GWR.
Attention, ce qui en théorie peut paraître assez tranché ne l’est souvent pas du tout en pratique. En effet, il y a bien des cas où l’on a du mal à savoir dans quel cadre on se situe exactement.
Pour réaliser une GWR sur R plusieurs packages reconnus existent. On
peut citer notamment le package spgwret le
package GWmodel. Nous choisirons d’utiliser ici le
package GWmodel.
5.1 Calcul matrice des distance
La première étape est de calculer la distance entre tous nos
observations. Pour ce faire nous utiliserons la fonction
gw.dist().
library(GWmodel)
# Le package GWmodel n'est pas compatible avec le format sf il a besoin d'un shape (contrairement à spgwr qui peut travailler avec un format sf)
# Pour construire la matrice de distances entre centroïdes des EPCI :
dm.calib <- gw.dist(dp.locat = coordinates(shape))5.2 Définition de la bande passante
La bande passante est une distance au-delà de laquelle le poids des
observations est considéré comme nul. Le calcul de cette de distance
très importante car la valeur de la bande passante pourra fortement
influencer notre modèle. La définition de la bande passante renvoit à
quel type de pondération nous souhaitons appliquer. Heureusement la
fonction bw.gwrva choisir pour nous le résultat
optimal…
Pour ce faire la fonction va se baser sur un critère statistique que l’utilisateur devra définir: le CV (validation croisée) ou le AIC (Critère d’information d’Akaike). Elle reposera aussi sur un type de noyau qu’il faudra également définir : Gaussien, Exponentielle, Bicarré, Tricube ou encore Boxcar. Enfin, il sera également nécessaire de savoir si ce noyau pourra être adaptatif ou fixe.
Voici quelques informations pour guider nos choix :
- Le critère CV a pour objectif de maximiser le pouvoir prédictif du modèle, le critère AIC va chercher un compromis entre le pouvoir prédictif du modèle et son degrée de complexité. En général, le critère AIC est privilégié.
- Avec un noyau fixe l’étendu du noyau est determiné par la distance au point d’intérêt et il est identique en tout point de l’espace. Un noyau fixe est adapté si la répartition des données est homogène dans l’espace, l’unité de la bande passante sera donc une distance. Avec un noyau adaptatif l’étendu du noyau est détermininé par le nombre de voisins. Il est donc plus adapté à une répartition non homogène, l’unité sera alors le nombre de voisins.
Concernant la forme des noyaux :
- Les noyaux gaussiens et exponentiels vont pondérer toutes les observations avec un poids qui tend vers zéro avec la distance au point estimé.
- Les noyaux Bisquare et Tricube (dont les forme est très proche) accordent également aux observations un poids décroissant avec la distance, mais par contre ce poids est nul au delà de la distance définit par la bande passante.
- Le noyau Box-Car traite un phénomène continu de façon discontinue.
Sachant que sur la forme du noyau, il est tout à fait possible de comparer deux pondérations et deux modèle de GWR.
# Définition de la bande passante (bandwidth en anglais) :
bw_g <- bw.gwr(data = shape,
approach = "AICc",
kernel = "gaussian",
adaptive = TRUE,
dMat = dm.calib,
formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)Adaptive bandwidth (number of nearest neighbours): 763 AICc value: 17959.04
Adaptive bandwidth (number of nearest neighbours): 480 AICc value: 17884.14
Adaptive bandwidth (number of nearest neighbours): 303 AICc value: 17798.07
Adaptive bandwidth (number of nearest neighbours): 196 AICc value: 17715.24
Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 17622.76
Adaptive bandwidth (number of nearest neighbours): 87 AICc value: 17526.89
Adaptive bandwidth (number of nearest neighbours): 59 AICc value: 17422.28
Adaptive bandwidth (number of nearest neighbours): 45 AICc value: 17353.68
Adaptive bandwidth (number of nearest neighbours): 33 AICc value: 17283.09
Adaptive bandwidth (number of nearest neighbours): 29 AICc value: 17261.83
Adaptive bandwidth (number of nearest neighbours): 23 AICc value: 17250.23
Adaptive bandwidth (number of nearest neighbours): 22 AICc value: 17247.23
Adaptive bandwidth (number of nearest neighbours): 19 AICc value: 17201.79
Adaptive bandwidth (number of nearest neighbours): 19 AICc value: 17201.79
bw_g[1] 19
5.3 Estimation du modèle
Une fois la bande passante définit on peut lancer l’estimation de notre modèle de GWR
mod.gwr_g <- gwr.robust(data = shape,
dMat = dm.calib,
bw = bw_g,
kernel = "gaussian",
filtered = FALSE, #un des problèmes de la GWR est de gérer des individus "aberrants" au niveau local. 2 méthodes ont été définit pour gérer cela.
# Méthode 1 (argument TRUE) on filtre en fonction des individus standardisés. L'objectif est de détecter les individus dont les résidus sont très élevés et de les exclure.
# Methode 2 (argument FALSE) on diminue le poids des observations aux résidus élevés.
adaptive = TRUE,
formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)Si on souhaite comparer deux modèles car nous avons un doute sur les paramètre c’est tout à fait possible. Par exemple ici nous souhaitons comparer deux formes de noyau :
bw_tri <- bw.gwr(data = shape,
approach = "AICc",
kernel = "tricube",
adaptive = TRUE,
dMat = dm.calib,
formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)Adaptive bandwidth (number of nearest neighbours): 763 AICc value: 17701.17
Adaptive bandwidth (number of nearest neighbours): 480 AICc value: 17586.94
Adaptive bandwidth (number of nearest neighbours): 303 AICc value: 17422.51
Adaptive bandwidth (number of nearest neighbours): 196 AICc value: 17294.75
Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 17254.1
Adaptive bandwidth (number of nearest neighbours): 87 AICc value: 17279.33
Adaptive bandwidth (number of nearest neighbours): 154 AICc value: 17259.05
Adaptive bandwidth (number of nearest neighbours): 112 AICc value: 17257.17
Adaptive bandwidth (number of nearest neighbours): 138 AICc value: 17254.8
Adaptive bandwidth (number of nearest neighbours): 122 AICc value: 17253.75
Adaptive bandwidth (number of nearest neighbours): 117 AICc value: 17255.04
Adaptive bandwidth (number of nearest neighbours): 123 AICc value: 17253.57
Adaptive bandwidth (number of nearest neighbours): 126 AICc value: 17254.25
Adaptive bandwidth (number of nearest neighbours): 123 AICc value: 17253.57
mod.gwr_tri <- gwr.robust(data = shape,
dMat = dm.calib,
bw = bw_tri,
kernel = "gaussian",
filtered = FALSE,
adaptive = TRUE,
formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)
Best_gwr <- cbind(
rbind(bw_g, bw_tri),
rbind(mod.gwr_g$GW.diagnostic$gw.R2,mod.gwr_tri$GW.diagnostic$gw.R2),
rbind(mod.gwr_g$GW.diagnostic$AIC,mod.gwr_tri$GW.diagnostic$AIC)) %>%
`colnames<-`(c("Nb Voisins","R2","AIC")) %>%
`rownames<-`(c("GAUSSIAN","TRICUBE"))
Best_gwr Nb Voisins R2 AIC
GAUSSIAN 19 0.9324328 16849.26
TRICUBE 123 0.8577370 17567.05
Le modèle avec une forme qui a été définit au format gaussien explique un meilleur \(R^2\) et le score d’\(AIC\) est plus faible. Ce modèle est donc plus performant et dans notre situation c’est plutôt sur ce modèle qu’il faut partir.
5.4 Interprétation des premiers résultats
Comme pour le modèle linéaire classique l’objet qui contient notre GWR est composé de plusieurs éléments. Pour obtenir nos résultats il suffit d’appeler l’objet.
# Pour voir les différent élément qui compose notre modèle de GWR
summary(mod.gwr_g) Length Class Mode
GW.arguments 11 -none- list
GW.diagnostic 8 -none- list
lm 14 lm list
SDF 1223 SpatialPolygonsDataFrame S4
timings 5 -none- list
this.call 13 -none- call
Ftests 0 -none- list
# Pour accéder aux résultat
mod.gwr_g ***********************************************************************
* Package GWmodel *
***********************************************************************
Program starts at: 2022-09-29 22:37:20
Call:
gwr.basic(formula = formula, data = data, bw = bw, kernel = kernel,
adaptive = adaptive, p = p, theta = theta, longlat = longlat,
dMat = dMat, F123.test = F123.test, cv = T, W.vect = W.vect)
Dependent (y) variable: prix_med
Independent variables: perc_log_vac perc_maison perc_tiny_log dens_pop med_niveau_vis part_log_suroccup part_agri_nb_emploi part_cadre_profintellec_nbemploi
Number of data points: 1223
***********************************************************************
* Results of Global Regression *
***********************************************************************
Call:
lm(formula = formula, data = data)
Residuals:
Min 1Q Median 3Q Max
-1566.8 -220.2 -27.2 174.0 3277.3
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1592.873 11.204 142.171 < 2e-16 ***
perc_log_vac -287.337 14.134 -20.330 < 2e-16 ***
perc_maison -128.255 20.041 -6.400 2.22e-10 ***
perc_tiny_log -261.556 27.934 -9.363 < 2e-16 ***
dens_pop 173.942 15.272 11.389 < 2e-16 ***
med_niveau_vis 288.017 13.860 20.781 < 2e-16 ***
part_log_suroccup 388.982 27.352 14.221 < 2e-16 ***
part_agri_nb_emploi -20.785 15.059 -1.380 0.168
part_cadre_profintellec_nbemploi -7.841 19.904 -0.394 0.694
---Significance stars
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 391.8 on 1214 degrees of freedom
Multiple R-squared: 0.7784
Adjusted R-squared: 0.7769
F-statistic: 532.9 on 8 and 1214 DF, p-value: < 2.2e-16
***Extra Diagnostic information
Residual sum of squares: 186354789
Sigma(hat): 390.6721
AIC: 18086.13
AICc: 18086.31
BIC: 16985.31
***********************************************************************
* Results of Geographically Weighted Regression *
***********************************************************************
*********************Model calibration information*********************
Kernel function: gaussian
Adaptive bandwidth: 19 (number of nearest neighbours)
Regression points: the same locations as observations are used.
Distance metric: A distance matrix is specified for this model calibration.
****************Summary of GWR coefficient estimates:******************
Min. 1st Qu. Median
Intercept 1109.13394 1309.99018 1516.79508
perc_log_vac -919.82407 -224.60028 -163.40316
perc_maison -898.16365 -258.64481 -110.20592
perc_tiny_log -1210.96882 -206.68578 -92.91874
dens_pop -411.20172 72.82448 217.03593
med_niveau_vis 81.29015 287.78992 358.32914
part_log_suroccup -543.68600 42.42839 142.48859
part_agri_nb_emploi -498.74706 -65.81443 -32.88823
part_cadre_profintellec_nbemploi -347.37935 -48.57329 0.58811
3rd Qu. Max.
Intercept 1696.07367 2330.78
perc_log_vac -113.83881 388.64
perc_maison 76.01079 896.95
perc_tiny_log 32.81544 773.30
dens_pop 268.35651 659.84
med_niveau_vis 413.36560 853.98
part_log_suroccup 247.33650 700.50
part_agri_nb_emploi -1.15602 120.33
part_cadre_profintellec_nbemploi 49.93194 388.50
************************Diagnostic information*************************
Number of data points: 1223
Effective number of parameters (2trace(S) - trace(S'S)): 320.246
Effective degrees of freedom (n-2trace(S) + trace(S'S)): 902.754
AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 17201.79
AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 16849.26
BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 17068.01
Residual sum of squares: 56808914
R-square value: 0.9324328
Adjusted R-square value: 0.9084372
***********************************************************************
Program stops at: 2022-09-29 22:38:08
Au niveau des résultats, on a d’abord un rappel complet du modèle linéaire classique. Puis vient ensuite les informations concernant notre GWR.
Ce que l’on peut relever immédiatement c’est que le \(R^2 ajusté\) de la GWR est nettement meilleur que celui de la régression linéaire multiple. On passe de \(77\%\) de variance expliqué à \(91\%\) avec la GWR. La seconde information qui nous intéresse particulièrement se sont les coefficient associés à nos VI. On voit immédiatement qu’ils ne sont pas présenté de la même manière que ceux de la régression linéaire. En effet, chaque VI va avoir des coefficent en fonction du minimum, maximum et des quartiles. Cela permet de rendre compte de la stationnarité de l’effet ou non. Dans notre cas on voit qu’il y a bien une variation et même dans certain cas une inversion des signes. Cela laisse supposer une non stationnarité des effets, c’est à dire qu’il y aurait un effet local qui ne suivrait pas l’effet global.
Par exemple, pour le pourcentage de logement vacant avec un coefficient global (modèle linéaire) de \(-287\), quand ce pourcentage augmente le prix médian baisse. En simplifiant le pourcentage baisse d’1% le prix median augmente de \(287€\) (si l’unité du prix médian est l’€). Dans le cas de la densité de population on a un coefficient global positif donc une relation positive. La densité augmente donc le prix médian augmente. Ici, au global la densité augmente de 1 le prix médian augmente de \(173€\)
Les résultats de la GWR illustrent comment les coefficient varient en fonction des unités spatiales, on est bien sur une approche locale. Gardons avec l’exemple de la densité de population. Dans les lieux où le prix médian est à son minimum le coefficient est de \(-411\); On a donc une relation négative. Dans ces epaces la densité augment de1 le prix médian baisse de \(411\). Ensuite on va constater une inversion du signe du coefficient. Ainsi dans les EPCI dans le dernier quartile où le prix médian du logement est le plus élevé (par ex Paris) le coefficient est positif. A son maximum une augmention d’1 unité entraine une augmentation du prix de \(659€\). On a donc très clairement un effet de la densité qui ne sera pas du tout le même en fonction du lieu. Ce qui est aussi très intéressant c’est qu’on peut étudier l’intervalle interquartile. Ainsi, toujours pour la densité, ça implique que pour 50% de nos unités spatiales (EPCI entre quartile 1 et 3) une augmentation d’une unité de la densité va impliquer une augmentation du prix médian entre \(72€\) et \(268€\).
Au travers de ces résultats on voit parfaitement comment une même variables peut avoir un effet différents voir opposer en fonction des unités de lieu.
La cartographie va être la meilleure manière de représenter les betas (coefficients) et les différents indicateurs fournis avec la GWR, cela nous permet de décrire plus finement et plus précisément les phénomènes observés.
L’ensemble des données est stocké dans le sous objet SDF de notre modèle. Il contient l’ensemble des informations du modèle associé à chaque donnée spatiale.
On peut le convertir en un dataframe pour le visualiser plus facilement. A l’origine il est au format “SpatialPointsDataFrame”.
# Pour visualiser ce fichier dans R
#View(mod.gwr_g$SDF@data)
#Pour voir à quoi il ressemble dans ce document
datatable(mod.gwr_g$SDF@data)# Pour voir les variables qui le constituent
names(mod.gwr_g$SDF@data) [1] "Intercept" "perc_log_vac"
[3] "perc_maison" "perc_tiny_log"
[5] "dens_pop" "med_niveau_vis"
[7] "part_log_suroccup" "part_agri_nb_emploi"
[9] "part_cadre_profintellec_nbemploi" "y"
[11] "yhat" "residual"
[13] "CV_Score" "Stud_residual"
[15] "Intercept_SE" "perc_log_vac_SE"
[17] "perc_maison_SE" "perc_tiny_log_SE"
[19] "dens_pop_SE" "med_niveau_vis_SE"
[21] "part_log_suroccup_SE" "part_agri_nb_emploi_SE"
[23] "part_cadre_profintellec_nbemploi_SE" "Intercept_TV"
[25] "perc_log_vac_TV" "perc_maison_TV"
[27] "perc_tiny_log_TV" "dens_pop_TV"
[29] "med_niveau_vis_TV" "part_log_suroccup_TV"
[31] "part_agri_nb_emploi_TV" "part_cadre_profintellec_nbemploi_TV"
[33] "E_weigts" "Local_R2"
# Intercept : c'est la constante c'est à dire prix médian de référence
# nom de la variable : estimation du coefficient, du beta associé à la VI en chaque point.
# y : les valeurs de la VD
# yhat : valeur de y prédite.
# residual, Stud_residual : résidu et résidu standardisé
# CV_score : score de validation croisée
# _SE : erreur standard de l’estimation du coefficient
# _TV : t-value de l’estimation du coefficient
# E_weight : poids des observations dans la régression robuste
# Local_R2 : R2 au niveau de chaque unité spatiale5.4.1 Etude des résidus
Commençons par une étude des résidus afin de vérifier que cette fois elles n’ont pas de structure apparente.
res_gwr <- mod.gwr_g$SDF$Stud_residual
data_immo$res_gwr <- res_gwr
# carto résidus v2
mf_map(x = data_immo,
var = "res_gwr",
type = "choro",
border = "#ebebeb",
lwd = 0.1,
breaks = quantile(data_immo$res_gwr, seq(0,1, by=1/11)),
pal = cols_v2,
leg_title = "Résidus GWR",
leg_val_rnd = 1)
mf_title("Résidus GWR")Il semble qu’il y ait bien une répartition aléatoire des résidus.
5.4.2 Etude des coefficients
Pour visualiser la non stationnarité des effets de nos VI la solution la plus efficace c’est la carte.
# On ajoute à data_immo les coefficients
data_immo$agri.coef=mod.gwr_g$SDF$part_agri_nb_emploi
data_immo$perc_maison.coef <- mod.gwr_g$SDF$perc_maison
data_immo$dens_pop.coef=mod.gwr_g$SDF$dens_pop
data_immo$med_vie.coef=mod.gwr_g$SDF$med_niveau_vis
data_immo$logvac.coef=mod.gwr_g$SDF$perc_log_vac
data_immo$tinylog.coef=mod.gwr_g$SDF$perc_tiny_log
data_immo$suroccup.coef=mod.gwr_g$SDF$part_log_suroccup
data_immo$cadre.coef=mod.gwr_g$SDF$part_cadre_profintellec_nbemploi
# Réaliser la collection des cartes
#par(mfrow = c(2, 4))
mf_map(x = data_immo, var = "agri.coef", type = "choro", pal= "Earth")
mf_title("Coefficients Agriculteurs")mf_map(x = data_immo, var = "perc_maison.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Maison")mf_map(x = data_immo, var = "dens_pop.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de dens pop")mf_map(x = data_immo, var = "med_vie.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Médianne niveau de vie")mf_map(x = data_immo, var = "logvac.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Logements vacants")mf_map(x = data_immo, var = "tinylog.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Petits logements")mf_map(x = data_immo, var = "suroccup.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de logement suroocupés")mf_map(x = data_immo, var = "cadre.coef", type = "choro", pal= "Earth")
mf_title("Coefficients de Cadre")Les cartes des betas vont ainsi illustrer la variation des effets en fonctions des entités spatiales. Dans notre cas qu’elles sont les EPCI où l’effet du coefficient est négatif et ceux où il est positif. Dans quel EPCI notre VI va entrainer une augmentatation du prix median dans quel autre au contraire une diminution. Sachant que dnas notre cas toutes les VI sont significative, elles ont donc toutes un effet qui varie localement
Par contre, ici on va étudier VI par VI, il peut être intéressant de savoir par EPCI quelle variable va le plus expliquer la variation de notre VD, laquelle a l’impact le plus important. C’est la carte des contributions max par EPCI. Pour la réaliser on va se baser sur le t-value.
data_immo$agri.t <- mod.gwr_g$SDF$part_agri_nb_emploi_TV
data_immo$maison.t <- mod.gwr_g$SDF$perc_maison_TV
data_immo$dens.t <- mod.gwr_g$SDF$dens_pop_TV
data_immo$medvie.t <- mod.gwr_g$SDF$med_niveau_vis_TV
data_immo$logvac.t <- mod.gwr_g$SDF$perc_log_vac_TV
data_immo$tinylog.t <- mod.gwr_g$SDF$perc_tiny_log_TV
data_immo$suroccup.t <- mod.gwr_g$SDF$part_log_suroccup_TV
data_immo$cadre.t <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi_TV
#Définir contrib max
df<- as.data.frame(data_immo)
# On passe les t-values en valuer absolues pour voir la plus grande contribution dans un sens sens ou ans l'autre
data_immo$contribmax<- colnames(df[, c(30:37)])[max.col(abs(df[, c(30:37)]),ties.method="first")]
par(mfrow = c(1, 1))
# Carte
mf_map(x = data_immo, var = "contribmax", type = "typo", pal= "Zissou 1")
mf_title("Carte des varibales contribuant le plus par epci")On note donc ici que dans le grand bassin parisien c’est la variable densité de population qui a l’effet le plus important sur le prix médian
On peut également cartographier les R2 locaux. Ce qui donnera une indication sur les zone où la variabilité est la mieux expliqué.
data_immo$r2local=mod.gwr_g$SDF$Local_R2
mf_map(x = data_immo, var = "r2local", type = "choro", pal= "Reds")
mf_title("R2 Locaux")A partir des t-value on peut aussi étudier la significativité des effets sur le territoire. On peut ainsi calculer et cartographier un indicateur qui représenterait le nombre de VI dont l’effet est significatif sur chaque unité spatiale. Cela donne une bonne idée de la complexité du phenomène sur un espace donnée (en effet sur un epci on peut avoir toutes les variables de signicatives, elle jouent donc sur cette espace toutes un rôles)et souligne l’importance d’avoir une carte par coefficient.
# Pour rappel si on a plus de 200 individus et le t-value > |1.96| on pourra considérer le coefficient comme significatif au seuil de 0.05 (95% chances que ce ne soit pas dû au hasard)
data_immo$nbsignif_t<- rowSums(abs(df[, c(30:37)]) > 1.96)
mf_map(x = data_immo, var = "nbsignif_t", type = "typo", pal= "Reds")
mf_title("Nombre des Betas significatif par EPCI (t-value)")Il se peut que cela soit plus intéressant d’utiliser les p-value, notamment si vous avez moins de 200 individus.
# Les p-value ne sont pas fournit dans le modèle de la GWR, on pourrait les calculer à partir de t-value et de l'erreur standard mais le package GWmodel propose une fonction pour les obttenir
pvalue<-gwr.t.adjust(mod.gwr_g)
# On ajoute les p-value à notre fichier
data_immo$agri.p <- pvalue$SDF$part_agri_nb_emploi_p
data_immo$maison.p <- pvalue$SDF$perc_maison_p
data_immo$dens.p <- pvalue$SDF$dens_pop_p
data_immo$medvie.p <- pvalue$SDF$med_niveau_vis_p
data_immo$logvac.p <- pvalue$SDF$perc_log_vac_p
data_immo$tinylog.p <- pvalue$SDF$perc_tiny_log_p
data_immo$suroccup.p <- pvalue$SDF$part_log_suroccup_p
data_immo$cadre.p <- pvalue$SDF$part_cadre_profintellec_nbemploi_p
df<- as.data.frame(data_immo)
data_immo$nbsignif_p<- rowSums(df[, c(41:48)] < 0.05)
mf_map(x = data_immo, var = "nbsignif_p", type = "typo", pal= "Reds")
mf_title("Nombre des Betas significatif par EPCI (p-value)")Il est aussi possible de réaliser une collection de carte des p-value (ou t-value) comme ce qui a été fait pour les coefficients. L’intérêt est de voir où l’effet de la VI est significatif et où il ne l’est pas.
# Ici nous représenterons les p-value avec un decoupage par classe de significativité et seulement les p-value de 2 VI
par(mfrow = c(1, 2))
# Par exemple les p-value des coefficients de la variable part de l'emploi agriculteur
data_immo<- data_immo %>% mutate(agri.p_fac = case_when(agri.p<= 0.002 ~ "[0;0.002[",
agri.p <= 0.01 ~ "[0.002;0.01[",
agri.p <= 0.05 ~ "[0.01;0.05[",
agri.p <= 0.1 ~ "[0.05;0.1[",
TRUE ~ "[0.1;1]")) %>%
mutate(agri.p_fac = factor(agri.p_fac,
levels = c("[0;0.002[", "[0.002;0.01[",
"[0.01;0.05[",
"[0.05;0.1[", "[0.1;1]")))
mypal2 <- mf_get_pal(n = 5, palette = "SunsetDark")
mf_map(x = data_immo,
var = "agri.p_fac",
type = "typo",
border = "grey3",
lwd = 0.1,
pal=mypal2,
leg_title = "Classe P-value")
mf_title("P-value du coefficient de la part d'emploi agriculteurs")
# Pour la densité de population
data_immo<- data_immo %>% mutate(dens.p_fac = case_when(dens.p<= 0.002 ~ "[0;0.002[",
dens.p <= 0.01 ~ "[0.002;0.01[",
dens.p <= 0.05 ~ "[0.01;0.05[",
dens.p <= 0.1 ~ "[0.05;0.1[",
TRUE ~ "[0.1;1]")) %>%
mutate(dens.p_fac = factor(dens.p_fac,
levels = c("[0;0.002[", "[0.002;0.01[",
"[0.01;0.05[",
"[0.05;0.1[", "[0.1;1]")))
mypal2 <- mf_get_pal(n = 5, palette = "SunsetDark")
mf_map(x = data_immo,
var = "dens.p_fac",
type = "typo",
border = "grey3",
lwd = 0.1,
pal=mypal2,
leg_title = "Classe P-value")
mf_title("P-value du coefficient de la densité de population")Ces cartes des p-value est particulièrement importante car nous donne les endroit où l’effet est signifiatif. En effet, on sait que la VI a effet global qui est significatif, qu’elle en plus une variabilité locale or localement elle n’est pas partout significative. Pour la part d’agriculteur dans l’emploi a un effet qui est significatif quasiment que dans le sud-est
Conclusion
Nous avons donc réalisé une analyse du prix médian de l’immobilier en France métropolitaine. LA GWR nous semble être une méthode particulièrement intéressante à la croisé de la statistique et de la géographie, et qui peut s’avérer très utile globalement en SHS pour essayer d’appréhender la complexité des phénomènes qui sont étudiés. En effet, nous avons ainsi vue la limite des statistique dites “classique” pour appréhender des phénomènes avec un ancrage spatiale, la GWR nous permets non seulement de dépasser cette limite de l’indépendance tout en s’intéressant à la variabilité des effets des variables en fonction de l’espace.
La GWR n’est bien sur pas la seul approche existante pour s’intéresser à l’aspect spatial de phénomène et variable sociale, il existe des modèle de régressions spatiales (SDEM, SDM, SAR…) mais également d’autres méthode comme l’analyse territoriale multiscalaire (MTA) qui peuvent également s’avérer extremement intéressante et riche. # Bibliographie {-}
Annexes
Info session
| setting | value |
|---|---|
| version | R version 4.2.1 (2022-06-23) |
| os | Ubuntu 18.04.6 LTS |
| system | x86_64, linux-gnu |
| ui | X11 |
| language | (EN) |
| collate | fr_FR.UTF-8 |
| ctype | fr_FR.UTF-8 |
| tz | Europe/Paris |
| date | 2022-09-29 |
| pandoc | 2.18 @ /usr/lib/rstudio/bin/quarto/bin/tools/ (via rmarkdown) |
| package | ondiskversion | source |
|---|---|---|
| car | 3.0.11 | CRAN (R 4.1.2) |
| carData | 3.0.5 | CRAN (R 4.2.0) |
| correlation | 0.8.0 | CRAN (R 4.1.3) |
| corrplot | 0.90 | CRAN (R 4.1.2) |
| digest | 0.6.29 | CRAN (R 4.2.0) |
| dplyr | 1.0.9 | CRAN (R 4.2.0) |
| DT | 0.19 | CRAN (R 4.1.1) |
| GGally | 2.1.2 | CRAN (R 4.2.0) |
| ggplot2 | 3.3.6 | CRAN (R 4.2.0) |
| gtsummary | 1.5.0 | CRAN (R 4.1.1) |
| GWmodel | 2.2.8 | CRAN (R 4.1.2) |
| here | 1.0.1 | CRAN (R 4.2.0) |
| mapsf | 0.3.0 | CRAN (R 4.1.1) |
| maptools | 1.1.4 | CRAN (R 4.2.0) |
| Matrix | 1.5.1 | CRAN (R 4.2.1) |
| plotly | 4.10.0 | CRAN (R 4.1.1) |
| RColorBrewer | 1.1.3 | CRAN (R 4.2.0) |
| Rcpp | 1.0.9 | CRAN (R 4.2.1) |
| rgeoda | 0.0.9 | CRAN (R 4.1.3) |
| rmarkdown | 2.14 | CRAN (R 4.2.0) |
| robustbase | 0.95.0 | CRAN (R 4.2.0) |
| rzine | 0.1.0 | gitlab (rzine/package@3909f2593c092aa9287d193859d7e4ca7b239a7a) |
| sf | 1.0.3 | CRAN (R 4.1.1) |
| sp | 1.4.7 | CRAN (R 4.2.0) |
| spatialreg | 1.2.3 | CRAN (R 4.2.0) |
| spData | 2.0.1 | CRAN (R 4.1.2) |
| spdep | 1.1.12 | CRAN (R 4.1.2) |
Citation
Auteur.e P, Auteur.e S (2021). “Titre de la fiche.” doi:10.48645/xxxxxx, https://doi.org/10.48645/xxxxxx,, https://rzine.fr/publication_rzine/xxxxxxx/.
BibTex :
@Misc{,
title = {Titre de la fiche},
subtitle = {Sous-Titre de la fiche},
author = {Premier Auteur.e and Second Auteur.e},
doi = {10.48645/xxxxxx},
url = {https://rzine.fr/publication_rzine/xxxxxxx/},
keywords = {FOS: Other social sciences},
language = {fr},
publisher = {FR2007 CIST},
year = {2021},
copyright = {Creative Commons Attribution Share Alike 4.0 International},
}